#!/usr/bin/env python
#EnvsNode.py

"""
gets the environments for the descendants of a node and counts branch lengths
for environmental combos

Revision History:

3-8-04: Cathy Lozupone: started writing code: revision from tree_comparisons.py

3-16-04: added unique fraction method

3-18-04 to 3-24-04: added and tested uniFracTreeEnv, _calc_rand_uniFracTreeEnv,
and _make_randomized                                                               

4-6-04: employed Fitch algorithm for parsimony distance analysis and
ancestral node estimation (assignInternalEnvsFitch)

4-8-04: wrote code for generating saturation curves

4-15-04: got rid of Anc_envs property and made Envs into a set object

4-23-04: parameterized the _p_value_tree function: set up
get_saturation_curve so that it would work for individual envs in a tree

9-17-04: added icr property to store the interconversion rate for a node

11-29-04: update makeEnvMatrix to have a RowOrder and ColOrder for use with
clustering

2-7-05: added code for jackknifing UPGMA clustered env trees

4-6-05: updated Envs property to be an array object for 1)faster 
processesing and 2)enhancements to the UniFrac algorithm to account
for abundance information

4-13-05: added functions for UniFrac abundance algorithms

2-19-06: removed copy_tree since it is slower than deepcopy, added
envTaxaCountArray because it is faster than envTaxaCount

2-23-06: changed phylogenticTestTreeString, uniFracTree, and uniFracWeightedTree to corrects 0.0 p-values

changed _p_value_pairwise so that it deals with 0.0 p-values and corrects for multiple comparison

changed evaluate_threshold_nodes so that it corrects for multiple comparisons an
d it returns p-value information for the significant nodes and all of the nodes.
 Note the in the past I returned just a list of nodes and the p value could be a
 ccessed from the node.p property. Now the node.p property is set to the uncorrec
 ted p_value and the corrected p_values are returned from the function. The sig_n
 odes that is returned is a list of (node, p-value) tuples as is the all_nodes

 changes _find_sig_nodes_a, find_sig_nodes_list_a,  _find_sig_nodes_b, and _find_
 sig_nodes_list_b so that they all correct for multiple comparisons and return re
 sults for sig_nodes and all_evaluated nodes.

3-8-06: edited evaluateBalancedUPGMA and associated methods so that they work by reassign environmental assignments to a constant tree rather than removing actual terminal nodes from the tree. This was done in order to 1) make it so that the tool can handle trees in which one sequence is assigned to multiple environments, 2) handle weighted jackknifing
*** note - this may behave sightly differently from the old tool because a
branch that leads from the root to a single node is accounted for whereas
in the past it was ignored***

Changed envTaxaDict so that it gets EnvsNames from a TerminalDescendant rather
than from the root. The allows it to be used on trees for which the Internal
Environments have not been assigned

3-8-06: changed envTaxaCountArray to envTaxaCountWeight since it can only do
Weighted counts

4-2-06: changed uniFracTreeEnv so that it reports <= with p-values when necessary

4-10-06: made changes so that evaluate_threshold_nodes performs randomizations
of all of the envs in the tree when evaluating a node. Involved changes to
_calc_random_parameter, _p_value_tree, _find_sig_nodes_list_b, and
get_corrected_ps. made it so that nodeG_fit sets the obs_exp_counts attribute
to the nodes.

4_18_06: made changes to Weighted UniFrac so that it normalizes distances

4-25-06 Micah Hamady: added setTerminalDescendantsDistance. performance related changes. reduced time for buildenvmatrix from 9 days to 4 hours for large dataset example with weighted, from a few days to about 10 minutes with unweighted. 

4-26-06: Cathy Lozupone: added make_tree_list, batch_pairwise_unifrac_values, \
        and uniFracP_uncertainty. These all are for calculating UniFrac values
        from many equally likely trees.
5-3-06: made changes to evaluate_threshold_nodes so that it takes Abundance
        into account when evaluating min_leafs. Also added some code for batch
        tree analysis

7-7-06: made changes to clusterEnvsUPGMA, envDistanceMatrix, clusterEnvsPCA, and
evaluateBalancedUPGMA so that normalization of Weighted UniFrac values is
parameterized

"""

from __future__ import division
from old_cogent.base.tree import PhyloNode
from old_cogent.base.dict2d import Dict2D
from old_cogent.parse.tree import DndParser
from old_cogent.parse.record_finder import TailedRecordFinder
from copy import deepcopy
from random import shuffle, randrange
from sets import Set, ImmutableSet
import re
from old_cogent.parse.record import RecordError
from old_cogent.maths.stats.test import G_fit_from_Dict2D, \
        chi_square_from_Dict2D, calc_contingency_expected, MonteCarloP, \
        t_paired
from old_cogent.cluster.metric_scaling import principal_coordinates_analysis, \
        output_pca
from old_cogent.cluster.UPGMA import inputs_from_dict2D, UPGMA_cluster
from Numeric import *
from resampler import JackknifeTreeSpeciesUPGMA, JackknifeEnvsUPGMA
from time import time
#==========================================================

BIG_NUM = 1e305

class EnvsNode(PhyloNode):
    def __init__(self, *args, **kwargs):
        """Overwrite the TreeNode init and create Envs, EnvsNames properties.

        self.Envs is a numeric array object to store the the environments 
        that are represented by each node. Usually, a 1 denotes the the Env
        is associated with that node and a 0 denotes that the Env is not. 
        Information on the relative contribution of each Env for each node
        can also be represented by storing numbers other than 0 and 1. 

        self.EnvsNames is an ordered list that stores the names of the envs
        that are associated with the positions in the self.Envs array

	self.p is a p-value associated with the node. (can be 
	calculated with different algorithms)

        """
        
        super(EnvsNode, self). __init__(*args, **kwargs)
        self.EnvsNames = None
        self.Envs = None
        self.p = None

#-------------------------------------------
#Data property assignments
    def setTerminalDescendantsDistance(self):
        """ Calculates distance from root to term desc

        sets property RootDistance

        will raise an exception if not called from root.
        """
        if self.Parent:
            raise ValueError, "Self not root node!"

        for node in self.TerminalDescendants:
            dist_from_root = node.BranchLength or 0
            p = node.Parent
            while p:
                dist_from_root += p.BranchLength or 0
                p = p.Parent
            node.RootDistance = dist_from_root

    def addEnvsToData(self):
        """appends the Envs in node.Envs to the Data property of terminal decs

        This is so that it is easier to evaluate unique branches when viewing a
        tree after running the p or unique fraction tests"""
        for node in self.TerminalDescendants:
            envs = ''
            envs_in_node = [node.EnvsNames[i] for i in nonzero(node.Envs)]
            for env in envs_in_node:
                if env:
                    envs = envs + '_' + env
            if not node.Data:
                node.Data = envs
            else:
                node.Data += envs
         
#-----------------------------------------------
#Env property assignments

    def clearTermEnvs(self):
        """Clears Envs property for all terminal descendants"""
        for term_dec in self.TerminalDescendants:
            term_dec.Envs = None

    def clearInternalEnvs(self):
        """Clears Envs property for all internal nodes and root"""
        for node in self.traverse():
            if node.Children:
                node.Envs = None
        
    def assignTerminalEnvs(self, envs, all_envs):
        """sets Envs property for terminal nodes using dict envs

        envs is a dictionary mapping the taxa name (Data property) to an
        array containing a 0 or 1 for the absence or presence of that node
        in each env associated with the tree

        all_envs is a sorted list of the envs in the tree that corresponds 
        to the arrrays in envs (made in make_env_dict_text)

        Looks at each terminal descendant in self. If the Data property of the
        node is in envs, sets the Envs property.

        """
        #set self.EnvsNames to a list of all of the envs in the tree
        self.EnvsNames = all_envs
        
        for node in self.TerminalDescendants:
            #set node.EnvsNames for each term dec and initialize a node.Envs
            #array
            node.EnvsNames = self.EnvsNames
            #populate the node.Envs array based on info in envs
            if node.Data in envs:
               node.Envs = envs[node.Data]
            else:
               node.Envs = zeros(len(self.EnvsNames),)
    
    def normalizeEnvs(self):
        """divides the env array by the number of sequences from that env

        must be run after assignTerminalEnvs and before assignInternalEnvs
        if UniFrac calculations are to be made taking abundance into account
        """
        #create a new array that is the sum of all of the Envs arrays
        array_sum = zeros(len(self.EnvsNames),Float)
        for node in self.TerminalDescendants:
            array_sum += node.Envs
        #replace 0s with 0.1 to avoid dividing by 0
        array_sum = choose(less(array_sum, 0.1), (array_sum, 0.1))
        #divide each node.Envs by the array_sum array
        for node in self.TerminalDescendants:
            node.Envs = node.Envs/array_sum
    
    def assignInternalEnvs(self, Weight=False):
        """recursive method that sets the Envs property for non-terminal nodes

        also sets the EnvsNames property
        
        The state of Weight determines whether or not abundance information
        will be taken into account.
        
        If Weight is False, Envs is set to the union of Envs from all of the 
        children. This means that a 1 will be set it the env is present
        amongst the children and a 0 will be set if it is not

        If Weight is True, Envs is set by summing the Envs arrays of all of the
        children.
        
        assignTerminalEnvs must be run first to set Envs and EnvsNames of
        the children.

        If Weight is True, normalizeEnvs must also be run first to account for 
        differences in sample number between envs
        """
        self.clearInternalEnvs()
        for node in self:
            node.assignInternalEnvs(Weight)
            #set the Envs property to an empty array of the appropriate 
            #size if it is not already set. Get info on the appropriate size
            #and on EnvsNames from the children
            if not self.Envs:
                self.EnvsNames = node.EnvsNames
                self.Envs = zeros(len(self.EnvsNames),Float)
            if not Weight:
                self.Envs = logical_or(self.Envs, node.Envs)
            else:
                self.Envs = self.Envs + node.Envs

    def assignInternalEnvsFitch(self, env1=None, env2=None):
        """sets the Anc_Env property for internal nodes and returns # changes

        Employs the Fitch algorithm which evaluates each internal node in a
        post-order tree traversal. First the intersection of the Envs of the
        children is made. If the intersection is not an empty set, the Envs
        property is set to the intersection. If it is an empty set, the Envs
        property is set to the union and a change is counted. The total changes
        is the parsimony evaluation of the total distance of the tree. Envs
        is a set object containing the possible ancestral states for the node"""
        
        self.clearInternalEnvs()
        #variable for storing the # of state changes in the tree
        pars_count = 0

        #setting traverse inputs to 'False, True' make a post-order traversal
        for node in self.traverse(False, True):
            if node:
                #make a set that is intersection of all children
                curr = node[0].Envs
                for child in node:
                    curr = logical_and(curr, child.Envs)
                if not curr:
                    for child in node:
                        curr = logical_or(curr, child.Envs)
                    pars_count += 1
                node.Envs = curr
                node.EnvsNames = node.Children[0].EnvsNames
        return pars_count

    def reassignAncestralEnvs(self, anc_list):
        """reassigns the Envs property of the Ancestors of a node

        this is used after removing a node from a tree. It is less time
        consuming than reassigning the Envs property for all of the nodes
        in the tree as was done earlier"""
        for anc in anc_list:
            #first clear the current Envs assignments by setting all values
            #in the array to 0
            anc.Envs[:] = 0
            #next reassign Envs as the union of the Envs in all of the Children
            #note that the Envs of all Children of a node must be assigned 
            #before assigning the node. This is accomplished because the 
            #Ancestor list is in order from the tips to the root
            for child in anc:
                anc.Envs = logical_or(child.Envs, anc.Envs)
        
    def hasEnv(self, env):
        """returns true if node has env in its Envs set"""
        if env not in self.EnvsNames:
            return 0
        else:
            return self.Envs[self.EnvsNames.index(env)]

    def envEmpty(self, x = None):
        """returns true if Envs set is empty

        x is empty variable used so that markTermDecs will work with
        both EnvCheck and hasEnv"""
        return not self.Envs

    def markTermDecs(self, f, env):
        """returns a list of indices for which f(TermDec, env) is True"""
        return[i for i, term_dec in enumerate(self.TerminalDescendants) if \
               f(term_dec, env)]

    def returnEnvTermDecs(self, env):
        """returns list of term_decs that are in an env"""
        return [i for i in self.TerminalDescendants if i.hasEnv(env)]

    def returnEmptyTermDecs(self, x = None):
        """returns a list of term_decs for which Envs is empty Set"""
        term_decs_i = self.markTermDecs(EnvsNode.envEmpty, x)
        return[self.TerminalDescendants[i] for i in term_decs_i]
        
    def envTaxaDict(self):
        """returns a dict mapping the envs to a list of the taxa in them"""
        #make dict and add an empty list for each env at the root
        env_taxa_dict = {}
        #make an array from concatenating Envs of all terminal descendants
        term_decs = self.TerminalDescendants
        term_decs_envs = [dec.Envs for dec in term_decs]
        term_envs = array(term_decs_envs)
        #populate env_taxa_dict based on array information
        for index, env in enumerate(term_decs[0].EnvsNames):
            env_array = term_envs[:,index]
            env_taxa_dict[env]  = [term_decs[i] for i in nonzero(env_array)]
        return env_taxa_dict

    def envTaxaCount(self):
        """returns a dict mapping the envs to the number of taxa in them"""
        env_taxa_dict = self.envTaxaDict()
        return dict([(key, len(env_taxa_dict[key])) for key in env_taxa_dict])

    def envTaxaCountWeight(self):
        """returns a dict mapping the envs to the # of taxa in them Weighted
        """

        self.clearInternalEnvs()
        self.assignInternalEnvs(Weight=True)
        return dict(zip(self.EnvsNames, self.Envs)) 
        
    def removeEnvTermNodes(self, env = None, number = 1):
        """removes number of terminal nodes from a specific env from tree
    
        If env is None (default), it randomly removes any terminal descendant
        """
        env_term_decs = self.getEnvTermDecs(env)
        shuffle(env_term_decs)
    
        for i in range(number):
            #select a terminal descendant to remove at random
            removed = env_term_decs.pop()
            rem_anc = removed.Ancestors
            self.removeNodeAndOnlyChildParents(removed)

            #update the Envs property for the tree
            self.reassignAncestralEnvs(rem_anc)
            
    def removeEnvAssignments(self, env, number, Weight):
        """Unassigns Terminal Descendants from a specific env"""
        env_term_decs = self.getEnvTermDecs(env)
        env_index = self.EnvsNames.index(env)
        if not Weight:
            shuffle(env_term_decs)
            for i in range(number):
                #select a terminal descendant to remove at random
                removed = env_term_decs.pop()
                removed.Envs[env_index] = 0.0
        else:
            env_term_decs_Abun = []
            for node in env_term_decs:
                Abun = int(node.Envs[env_index])
                add = [node] * Abun
                env_term_decs_Abun.extend(add)
            shuffle(env_term_decs_Abun)
            for i in range(number):
                removed = env_term_decs_Abun.pop()
                Abun = node.Envs[env_index]
                removed.Envs[env_index] -= 1
    
    def removeAllEnvTermDecs(self, envs_to_remove):
        """removes taxa that have a particular env in Envs from the tree

        envs_to_remove is a list of envs for which taxa will be removed"""

        nodes_to_remove = []
        env_taxa_dict = self.envTaxaDict()
        for env in envs_to_remove:
            if env in env_taxa_dict:
                nodes_to_remove += env_taxa_dict[env]
        for node in nodes_to_remove:
            self.removeNodeAndOnlyChildParents(node)

    def keepEnvTermDecs(self, envs_to_keep):
        """removes all taxa that are not in any of the envs in envs_to_keep

        envs_to_keep is a list of envs"""

        envs_to_keep = Set(envs_to_keep)
        for node in self.TerminalDescendants:
            envs = [node.EnvsNames[i] for i in nonzero(node.Envs)]
            envs = Set(envs)
            if not envs_to_keep.intersection(envs):
                self.removeNodeAndOnlyChildParents(node)
                
    def balanceTreeEnvs(self, num_to_keep, Weight):
        """Balances the number of Envs in the tree by keeping num_to_keep
        
        num_to_keep: the number of individuals from each env to keep
        Weight: Whether to remove in an unweighted or weighted fasion
            if True: removes the specified number of total sequences
            if False: removes the specified number of unique sequences (so
                if a Node was observed 5 times, it removes all 5 observations
        """
        if Weight:
            counts = self.envTaxaCountWeight()
        else:
            counts = self.envTaxaCount()
        for env in counts:
            if counts[env] > num_to_keep:
                num_to_remove = counts[env] - num_to_keep
                self.removeEnvAssignments(env, int(num_to_remove), Weight)
    
    def pruneTreeEnvs(self, env1, env2):
        """prunes self to only have taxa in the environment in env1 and env2"""
        env_taxa_dict = self.envTaxaDict()
        all_nodes = self.TerminalDescendants
        #get a set of nodes to remove
        env1_nodes = env_taxa_dict[env1]
        env2_nodes = env_taxa_dict[env2]
        remove_nodes = []
        for node in all_nodes:
            if node not in env1_nodes and node not in env2_nodes:
                remove_nodes.append(node)
        for node in remove_nodes:
            self.removeNodeAndOnlyChildParents(node)

    def getEnvTermDecs(self, env):
        """returns a list of term decs for env or all term decs if env is None
        """
        term_decs = self.TerminalDescendants
        #set env_term_decs to the full list of terminal descendants if env is
        #None otherwise set to the terminal descendants for that env only 
        if env == None:
            env_term_decs = term_decs
        else:
            env_term_decs = self.returnEnvTermDecs(env)

        return env_term_decs

    def getEnvPairs(self):
        """yields all possible pair combinations of environments in the tree"""
        #get a list of all of the envs in the tree (Envs property of root)
        all_envs = [self.EnvsNames[i] for i in nonzero(self.Envs)]
        num_envs = len(all_envs)

        #yields each env combo
        for index, env in enumerate(all_envs):
            env1 = env
            index2 = index + 1
            while index2 < num_envs:
                env2 = all_envs[index2]
                index2 += 1
                yield env1, env2
        
    def findSpecialistNodes(self, Envs, frac=1.0):
        """returns a list of nodes that have frac of desc from the Envs in Envs

        Envs is a list of Environments that the node must be a specialist for.
        Function only returns the deepest nodes with this property. It will not
        return the childen since all of the children must also be specialists
        and can thus be accessed through the parent

        frac is the fraction of descendants that must be from Envs for the
        node to be considered a specialist
        
        """
        self.clearInternalEnvs()
        #assign the internal Envs with weighted so that each number in array
        #represents the number of decendants from that Env (not presence/abs)
        self.assignInternalEnvs(Weight=True)
        specialists = []
        for node in self.traverse():
            if node.Parent:
                assigned_dec = sum(node.Envs)
                env_dec = sum([node.Envs[node.EnvsNames.index(env)] \
                        for env in Envs])
                if not assigned_dec:
                    spec_frac = 0.0
                else:
                    spec_frac = env_dec/assigned_dec
                if spec_frac >= frac and node.Parent not in specialists:
                    specialists.append(node)
                    
#                specialist_Envs = zeros(len(node.Envs), Float)
#                for env in Envs:
#                    specialist_Envs[node.EnvsNames.index(env)] = 1
#                if list(node.Envs) == list(specialist_Envs) and \
#                        list(node.Parent.Envs) != list(specialist_Envs):
#                    specialists.append(node)
        return specialists
                
    def expand_OTUs(self):
        """expands Terminal Descendants that represent more than 1 clone

        whether a sequence that was represented more than once is in the
        tree more than once affects randomization and jackknifing. This
        function creates a new node for each time a sequence was observed 
        and saves it as a child of the first node with a branch length of 0
        
        Tree should have environments assigned but not normalized"""
    
        for node in self.TerminalDescendants:
            if sum(node.Envs) > 1:
                num_envs = len(node.Envs)
                name_index = 0
                for i, env_count in enumerate(node.Envs):
                    if env_count > 0:
                        for count in range(int(env_count)):
                            new_node = EnvsNode()
                            new_node.Data = node.Data + '_s' + str(name_index)
                            name_index += 1
                            new_node.BranchLength = 0.0
                            new_node.Envs = zeros(num_envs,)
                            new_node.Envs[i] = 1.0
                            new_node.EnvsNames = node.EnvsNames
                            new_node.Parent = node
            
#--------------------------------------
#branch length evaluation for UniFrac

    def envBranchLengthsTree(self, results=None):
        """returns the branch length for every possible combo of envs
        
        Envs property should be previously defined with assignInternalEnvs"""
        #make initial dict to populate
        results = {}
        #make key for each individual env and set to 0 (in case nothing unique)
        envs = [self.EnvsNames[i] for i in nonzero(self.Envs)]
        for env in envs:
            results[ImmutableSet([env])] = 0
        #do a post-order tree traversal
        for node in self.traverse(self_before = False, self_after = True):
            if node.Parent is not None: #not at root
                #key for results dict is node.Envs as an immutable set
                node_envs = [node.EnvsNames[i] for i in nonzero(node.Envs)]
                envs_key = ImmutableSet(node_envs)
                #only count bl if the envs have been defined
                if envs_key and envs_key not in results:
                    results[envs_key] = 0
                if envs_key:
                    #add the branch length to appropriate key
                    results[envs_key] += node.BranchLength
        return results

    def envBranchLengthsPairwise(self):
        """calculates the unique and shared branchlengths for pairs of envs

        comboLengths maps the envs to the branchlength that includes
        them"""
        totalLengths = self.envBranchLengthsTree()
        #get each pair
        for env1, env2 in self.getEnvPairs():
            yield self.envBranchLengthsPair(env1, env2, totalLengths)

    def envBranchLengthsPair(self, env1, env2, totalLengths):
        """calculates the unique and shared branchlengths for a single pair of envs

        comboLengths maps the envs to the branchlength that includes
        them"""
        #initialize dict to store values
        comboLengths = {ImmutableSet([env1]): 0,
                        ImmutableSet([env2]): 0,
                        ImmutableSet([env1, env2]): 0
                        }
        #for each key see if is unique to env1, env2 or shared and add
        #bl to approprate category
        for key in totalLengths:
            #only in env1
            if env1 in key and env2 not in key:
                comboLengths[ImmutableSet([env1])] += totalLengths[key]
            #only in env2
            elif env1 not in key and env2 in key:
                comboLengths[ImmutableSet([env2])] += totalLengths[key]
            #in both env1 and env2
            elif env1 in key and env2 in key:
                comboLengths[ImmutableSet([env1, env2])] += \
                    totalLengths[key]
        return comboLengths

    def envUniqueFraction(self):
        """Maps each env to the fraction of the tree unique to that env."""
        self.clearInternalEnvs()
        self.assignInternalEnvs()
        lengths = self.envBranchLengthsTree()
        env_fractions = self._env_unique_fraction(lengths)
        return env_fractions

    def envUniqueFractionPairwise(self):
        """gets unique fractions for pairwise env comparisons"""
        self.clearInternalEnvs()
        self.assignInternalEnvs()
        for comboLengths in self.envBranchLengthsPairwise():
            env_fractions = self._env_unique_fraction(comboLengths)
            yield env_fractions

    def _env_unique_fraction(self, lengths_info):
        """calcs the fraction of the tree unique to each env using lengths_info

        lengths_info is generated with envBranchLengthsTree() if looking at the 
        whole tree or with envBranchLengthsPairwise() if making pairwise 
        comparisons of environments
        """
        env_fractions = {}
        #make sure that totalBranchLength is not 0 since dividing by it
        totalBranchLength = sum(lengths_info.values())
        if totalBranchLength:
            for key in lengths_info:
                if len(key) == 1:
                    env_fractions[key] = \
                                lengths_info[key] / totalBranchLength
        else:
            for key in lengths_info:
                if len(key) == 1:
                    env_fractions[key] = 0.0 
 
        return env_fractions
                    
    def totalUnique(self, env1=None, env2=None):
        """returns the total unique fraction in the tree
        
        env1 and env2 are placeholders so that can be used in general functions
        
        """
        env_fracs = self.envUniqueFraction()
        total_unique = sum(env_fracs.values())
        return total_unique

    def pairwiseTotalUnique(self):
        """gets fraction of tree unique to any env for all env combos"""
        for env_fracs in self.envUniqueFractionPairwise():
            envs = env_fracs.keys()
            total_unique = sum(env_fracs.values())
            #yields both the total fraction unique and the envs it is for
            yield list(envs[0])[0], list(envs[1])[0], total_unique

    def totalUniqueWeight(self, env1, env2):
        """calculates a distance between the two envs weighted for abundance"""
        self.normalizeEnvs()
        self.clearInternalEnvs()
        self.assignInternalEnvs(Weight=True)
        self.setTerminalDescendantsDistance()
        return self._total_unique_weight(env1, env2)
    
    def pairwiseTotalUniqueWeight(self, norm=False):
        """calculates pairwise distances between envs weighted for abundance

        calculated the distance between envs A and B as:
        d(A,B) = sum of abs(Ai - Bi) * branchlength for all nodes i

        assignTerminalEnvs, normalizeEnvs, and assignInternalEnvs(Weight = True)
        must be run first
        """
        self.normalizeEnvs()
        self.clearInternalEnvs()
        self.assignInternalEnvs(Weight=True)
        self.setTerminalDescendantsDistance()
        for env1, env2 in self.getEnvPairs():
            Distance = self._total_unique_weight(env1, env2, norm)
            yield env1, env2, Distance

    def _total_unique_weight(self, env1, env2, norm=False):
        """calculates the UniFrac distance between env1 and env2 weighted

        norm: if True, normalizes the distance for comparison between trees 
        by dividing them by the weighted average distance of each terminal 
        descendant from the root
        """
        #get the index for the environments
        st = time()
        env1i = self.EnvsNames.index(env1)
        env2i = self.EnvsNames.index(env2)
        Distance = 0
        st = time()
        for node in self.traverse():
            d = abs(node.Envs[env1i] - node.Envs[env2i])
            if node.BranchLength:
                Distance += d * node.BranchLength

        st = time()
        if norm:
                #normalize distance value
                avg_dist_from_root = 0
                for node in self.TerminalDescendants:
                #dist_from_root = self.distance(node)
                        dist_from_root = node.RootDistance
                        norm_dist_from_root = dist_from_root * \
                              (node.Envs[env1i] + node.Envs[env2i])
                        avg_dist_from_root += norm_dist_from_root
                Distance = Distance/avg_dist_from_root
        return Distance

    def totalEnvLengths(self):
        """returns the total branch lengths that include each environment

        total_lengths maps the env to the total branchlength that includes
        them"""
   
        total_lengths = {}
        lengths_info = self.envBranchLengthsTree()
        envs = [self.EnvsNames[i] for i in nonzero(self.Envs)]
        for env in envs:
            total_lengths[env] = 0
            for key in lengths_info:
                keys = list(key)
                if keys.count(env) >= 1:
                    total_lengths[env] += lengths_info[key]
        return total_lengths
            
    def makeEnvMatrix(self, default=None, Weight=False, norm=False):
        """makes a dict2D with totalUnique for all pairwise env comparisons
        
        also sets RowOrder property of the dict2D and Pads with default
        input
        
        Weight: if False, uses pairwiseTotalUnique which does not account for
        abundance information. If true, uses pairwiseTotalUniqueWeight which 
        does account for abundance information

        norm: if True, normalizes Weighted UniFrac values so that are a value
        between 0 and 1. If False, does not.
        """
        
        matrix = []
        if not Weight:
            for env1, env2, total_unique in self.pairwiseTotalUnique():
                matrix.append((env1, env2, total_unique))
                matrix.append((env2, env1, total_unique))
        else:
            for env1, env2, distance in self.pairwiseTotalUniqueWeight(norm):
                matrix.append((env1, env2, distance))
                matrix.append((env2, env1, distance))
        envs = [self.EnvsNames[i] for i in nonzero(self.Envs)]
	return Dict2D(matrix, RowOrder = envs, ColOrder = \
                envs, Pad=True, Default=default)

#===================================================
#Chi-square/ G test functions

    def makeEnvContingencyMatrix(self, node, Weight=True):
        """makes a Matrix representing an env contingency table for a node

        The Matrix has one row for each environment and 2 columns:
        The first column is the values for the node
        The second column is the values for the tree with the node removed
        """
        #check that not at the root
        if node.Parent:
            parent = node.Parent
            node.Parent = None
            if not Weight:
                result = {node.Data: node.envTaxaCount(),
                      'rest_of_tree': self.envTaxaCount()}
            else:
                result = {node.Data: node.envTaxaCountWeight(),
                      'rest_of_tree': self.envTaxaCountWeight()}

            result = Dict2D(result, Default=0,Pad=True)
            parent.append(node)
            return result

    def nodeChiSquare(self, node):
        """calculates the Chi Square p-value for a node"""

        obs_matrix = self.makeEnvContingencyMatrix(node)
        data_matrix = calc_contingency_expected(obs_matrix)
        test, node.p = chi_square_from_Dict2D(data_matrix)

    def nodeG_fit(self, node, Weight=True):
        """calculated the G_fit p-value for a node"""
        obs_matrix = self.makeEnvContingencyMatrix(node, Weight)
        data_matrix = calc_contingency_expected(obs_matrix)
        setattr(node, 'obs_exp_counts', data_matrix[node.Data])
        g_val, node.p = G_fit_from_Dict2D(data_matrix)
       
#=================================
#Calculating a P-value for the whole tree

    def uniFracTree(self, pop_size, node=None):
        """Computes one p- value for the tree with the UniFrac metric.
        
        if node is set calculates the p-value for the node but randomizes
        within self"""
        self._p_value_tree(EnvsNode.totalUnique, pop_size, node=node)
        raw_p = self.p
        corrected_p = str(raw_p)

        if self.p == 0.0:
            corrected_p = '<=' + str(1/float(pop_size))
            return "Corrected P Value: %s, Raw P Value: %f"  \
                            % (corrected_p, raw_p)

        return "P Value: %f"  % raw_p

    def uniFracWeightedTree(self, pop_size, node=None):
        """Computes tree p-value with weighted  UniFrac.
        
        the tree can only have sequences assigned to 2 different envs
        
        if node is set calculates the p-value for the node but randomizes
        within self"""
        #assign Internal Envs so that can determine which envs have sequences
        self.clearInternalEnvs()
        self.assignInternalEnvs()
        all_envs = [self.EnvsNames[i] for i in nonzero(self.Envs)]
        if len(all_envs) != 2:
            return 'NA'
        else:
            self._p_value_tree(EnvsNode.totalUniqueWeight, pop_size, \
                    tail='high', env1=all_envs[0], env2=all_envs[1], node=node)
            if self.p == 0.0:
                return '<=' + str(1/float(pop_size))
            else:
                return str(self.p)
        
    def phylogeneticTestTree(self, pop_size, node=None):
        """Computes one p - value for the tree with the Phylogenetic (P) test.

        If node is set calculates the p-value for the node but randomizes 
        within self.
        """
        self._p_value_tree(EnvsNode.assignInternalEnvsFitch, pop_size, 'low',
                node=node)
        raw_p = self.p
        corrected_p = str(self.p)

        if self.p == 0.0:
            corrected_p = '<=' + str(1/float(pop_size))
            return "Corrected P Value: %s, Raw P Value: %f"  \
                            % (corrected_p, raw_p)

        return "P Value: %f"  % raw_p

    def _p_value_tree(self, algorithm, pop_size, tail = 'high', env1=None, \
            env2=None, node=None):
        """Calculates a single p-value for the whole tree and assigns to self.p.

        The returned p-value indicates the collective probability that all 
        of the environments represented are not sampled from the same 
        distribution

        algorithm choices include: EnvsNode.totalUnique for UniFrac;
                        EnvsNode.totalUniqueWeight for weighted UniFrac;
                        EnvsNode.assignInternalEnvsFitch for Phylogenetic test

        pop_size is the number of random trees used to generate the random dist

        tail indicates which side of the distribution to look at:
            low = look for smaller values than expected by chance
            high = look for larger values than expected by chance

        if node is set, calculates the p_value for that node but does the 
            randomizations for the whole tree

        """
        #calculate value for true tree
        if not node or not node.Parent:
            value = algorithm(self, env1=env1, env2=env2)
        else:
            parent = node.Parent
            node.Parent = None
            value = algorithm(node, env1, env2)
            node.Parent = parent
            
        #calculate value for randomized trees
        random_tree = deepcopy(self)
        rand_values = random_tree._calc_rand_parameter(algorithm, pop_size, \
                env1=env1, env2=env2, node=node)

        if not node:
            self.p = MonteCarloP(value, rand_values, tail)
        else:
            node.p = MonteCarloP(value, rand_values, tail)
            

    def _calc_rand_parameter(self, algorithm, pop_size, env1=None, env2=None, \
            node=None):
        """calcs a value for the specified number of random trees

        returns a sorted list of values for random trees that are calculated
        based on the algorithm given

        Algorithm choices = self.totalUnique or self.assignInternalEnvsFitch
        or self.totalUniqueWeight
    
        method listed in algorithm should return a single value for each tree
        
        if node is set to a node other than the root, randomizes for the
        whole tree but returns the value for the node
        """

        values = []
        for i in range(pop_size):
            self._make_randomized()
            if not node or not node.Parent:
                values.append(algorithm(self, env1, env2))
            else:
                parent = node.Parent
                node.Parent = None
                values.append(algorithm(node, env1, env2))
                node.Parent = parent
        values.sort()
        return values

    def _make_randomized(self):
        """randomizes the terminal descendents info of self"""
        term_decs = self.TerminalDescendants
        #make a shuffled list of the envs for each terminal descendant
        envs = [node.Envs for node in term_decs]
        shuffle(envs)
        #reassign Envs for terminal nodes to Envs in shuffled list:
        self.clearTermEnvs()
        for i, node in enumerate(term_decs):
            node.Envs = envs[i]
    
#======================================
#calculating pairwise p values
    
    def uniFracPairwise(self, pop_size, Weight):
        """Uses the UniFrac metric to calculate pairwise P vals
        
        If Weight is True: uses weighted analysis 
        """
        if Weight:
            return self._p_value_pairwise(EnvsNode.totalUniqueWeight, \
                    pop_size, 'high', prune=True)
        else:
            return self._p_value_pairwise(EnvsNode.totalUnique, pop_size, \
                    'high', prune=True)
    
    def phylogeneticTestPairwise(self, pop_size):
        """Uses the Phylogenetic test to calculate pairwise P values"""
        return self._p_value_pairwise(EnvsNode.assignInternalEnvsFitch,\
                pop_size, 'low', prune=True)
            
    def _p_value_pairwise(self, algorithm, pop_size, tail ='high', prune=True):
        """returns a matrix with p-values for each pair of envs in a tree

        algorithm: the algorithm use to compute the p_value
            EnvsNode.totalUnique for UniFrac
            EnvsNode.totalUniqueWeight for weighted UniFrac
            EnvsNode.assignInternalEnvsFitch for the Phylogenetic test
            
        pop_size is the number of random trees used to generate the random dist

        tail indicates which side of the distribution to look at:
            low = look for smaller values than expected by chance
            high = look for larger values than expected by chance

        prune determines whether to remove the other envs from the tree first
            or to calc the p-values directly. The tree must be pruned for
            the phylogenetic test but not for UniFrac
        """

        matrix = []
        self.assignInternalEnvs()
        #get the number of comparisons being made
        num_envs = len([self.EnvsNames[i] for i in nonzero(self.Envs)])
        num_comp = sum([i for i in range(1,num_envs)])
        #get p-values for all pair combinations
        for env1, env2 in self.getEnvPairs():
            env_tree = deepcopy(self)
            if prune:
                env_tree.pruneTreeEnvs(env1, env2)
            env_tree._p_value_tree(algorithm, pop_size, tail, env1=env1, env2=env2)
            p_val = env_tree.p
            if p_val == 0.0:
                p_val = (1/float(pop_size)) * num_comp
                if p_val > 1.0:
                    p_val = 1.0
                p_val = '<=' + str(p_val)
            else:
                p_val = p_val * num_comp
                if p_val > 1.0:
                    p_val = 1.0
                p_val = str(p_val)

            matrix.append([env1, env2, p_val])
            matrix.append([env2, env1, p_val])
        return Dict2D(matrix, Pad=True)

#=================================
#calculating a separate p-value for each environment

    def uniFracTreeEnv(self, pop_size):
        """returns a dictionary of p-values: one for each env

        p-values indicate the probability that
        each environment is sampled from a unique distribution (is different
        from the collective others). The probability is the fraction of the
        unique branch lengths calculated from randomized trees that are greater
        than the observed unique branch lengths"""
        p_dict = {} #dictionary mapping environment to corrected p-value
        raw_p_dict = {} #dictionary mapping environment to raw p-value

        #calculate env fractions for true tree
        tree_fracs = self.envUniqueFraction()

        #calc env fractions for randomized trees
        #(returned as dict of sorted lists)
        random_tree = deepcopy(self)
        random_fracs = random_tree._calc_rand_uniFracTreeEnv(pop_size)

        #find the index for each env at which the tree unique fractions
        #is greater than the randomized unique fractions
        for key, real_frac in tree_fracs.items():
            p_val = MonteCarloP(real_frac, random_fracs[key])
            raw_p_dict[list(key)[0]] = str(p_val)
            if p_val == 0.0:
                p_val = '<=' + str((1/float(pop_size)))
            p_dict[list(key)[0]] = str(p_val)
        return p_dict, raw_p_dict

    def _calc_rand_uniFracTreeEnv(self, pop_size):
        """calcs the unique fracs for the specified number of randomized trees.

        returns a dict with sorted lists of branch lengths that are unique to
        each environment

        used with uniFracTreeEnv 
        """
        len_info = {}
        for i in range(pop_size):
            self._make_randomized()
            env_fracs = self.envUniqueFraction()
            for key in env_fracs:
                if key not in len_info:
                    len_info[key] = []
                    len_info[key].append(env_fracs[key])
                else:
                    len_info[key].append(env_fracs[key])
        for val in len_info.values():
            val.sort()
        return len_info

#======================================
    
    def clusterEnvsUPGMA(self, Weight=False, norm=False):
        """returns a tree clustering the environments in self with UPGMA.
        
        Weight: Determines whether to use normal(False) or Weighted(True)
                UniFrac
               
        norm: Determine whether to normalize the raw weighted UniFrac values
                (True) or not (False)"""

        #make a distance matrix for envs and pad with a large number
        env_matrix = self.makeEnvMatrix(BIG_NUM, Weight, norm)
        #turn dict2D into Numeric array and make list of phylonode objects
        matrix_a, node_order = inputs_from_dict2D(env_matrix)
        return UPGMA_cluster(matrix_a, node_order,BIG_NUM)

    def clusterEnvsPCA(self, Weight=False, norm=False):
        """Principal coordinates analysis of envs in self using UniFrac dist.

        returns tab delimited text that can be analyzed in Excel        
        """
        dis_mat = self.makeEnvMatrix(0.0, Weight, norm)
        matrix, node_order = inputs_from_dict2D(dis_mat)
        names = [node.Data for node in node_order]
        point_matrix, eigvals = principal_coordinates_analysis(matrix)
        return output_pca(point_matrix, eigvals, names)

    def jackknifeEnvsUPGMA(self, num_reps=100, num_envs=1, Weight=False):
        """measures robustness of nodes in UPGMA_cluster to removed Envs
        
        produces jackknife values for all nodes in a UPGMA_cluster by
        removing envs randomly from the distance matrix 
        and checking to see if the node is still supported

        self is a species tree (not an environment tree!)

        num_reps is the number of times to remove the envs when computing the
        jackknife values

        num_envs is the number of envs to remove

        Weight : if True, takes Abundance info into account, if False does not

        """
        jack = JackknifeEnvsUPGMA(self, num_reps, num_envs)
        return jack()

    def jackknifeSpeciesCurveUPGMA(self, env, stride=2, num_reps=100, \
            verbose=False, Weight=False):
        """measures robustness of the env-tree to # of species in env
        
        For the specified number of replicates (num_reps),
        removes the specified number (stride) of species from the 
        specified environment (env), reclusters the environments in the
        tree with UPGMA, and reports the % of the time that the new tree
        did not change for each number of removed environments.

        if verbose == True, prints the results to monitor progress"""
        
        tree2 = deepcopy(self)
        envTermDecs = tree2.getEnvTermDecs(env)
        num_taxa = len(envTermDecs)
        result = {}
        #remove species from the tree until there are less or equal species
        #in the tree than would would be removed in the next iteration
        while num_taxa > stride:
            #label each result entry based on the # of species in the tree
            num_taxa = num_taxa - stride
            jackknifer = JackknifeTreeSpeciesUPGMA(tree2, num_reps, \
                    stride, env=env, Weight=Weight)
            result[str(num_taxa)] = jackknifer.__call__()
            #remove env nodes in a copy of the original tree
            tree2.removeEnvTermNodes(env, stride)
            if verbose:
                print str(num_taxa), ':\t', result[str(num_taxa)]
        return result
                
    def jackknifeSpecies(self, num_removed, env, num_reps=100):
        """measures robustness of nodes to removed species from an env
        
        produces jackknife values for the nodes in an env cluster tree
        such as one made with clusterEnvsUPGMA
        
        num_removed is the number of species removed from the tree from the
            specified environment
        
        num_reps is the number of replicates
        
        env is the environment that the species will be removed from
        
        returns a dictionary containing the % of the time that a node is 
        supported when the species are removed"""

        result = {}
        #make a bit representation of the true UPGMA tree
        cluster = self.clusterEnvsUPGMA()
        cluster.nameUnnamedNodes()
        clusterArray, node_list = cluster.makeTreeArray()
        #initialize result array    
        for node in node_list:
            result[node.Data] = []
            
        #calculate jacknife values for all nodes 
        for i in range(num_reps):
            #remove the specified number of nodes from the specified 
            #environment in a copy of the original tree
            tree2 = deepcopy(self)
            tree2.removeEnvTermNodes(env, num_removed)
            cluster2 = tree2.clusterEnvsUPGMA()
            clusterArray2, node_list2 = cluster2.makeTreeArray()
        #compare array1 and array2. append a 1 to the result if a node in
        #the original tree is in the modified tree and a 0 if it is not
            for i_r, row in enumerate(clusterArray):
                    in_array = 0
                    for row2 in clusterArray2:
                        if list(row) == list(row2):
                            in_array = 1
                    result[node_list[i_r].Data].append(in_array)
        return result

    def evaluateBalancedUPGMA(self, num_reps, num_to_keep, Weight=False, \
            norm=False, verbose=False):
        """evaluate Balanced determines supported nodes in balanced env treesi
        
        num_reps: number of jackknife replicates to perform
        num_to_keep: The number of sequences to keep during each jackknife rep
        Weight: if False, makes tree using unweighted UniFrac; if True, uses
                weighted UniFrac
        norm: Whether or not to normalize the weighted UniFrac values by
                dividing by the average distance of seqs from the root
        """
        self.clearInternalEnvs()
        self.assignInternalEnvs(Weight)
        result = {}
        #record original env assignments for the Terminal Descendants
        all_envs = self.TerminalDescendants[0].EnvsNames
        env_dict = {}
        for node in self.TerminalDescendants:
            env_dict[node.Data] = node.Envs
        #do jackknifing
        for i in range(num_reps):
            jackknife_tree = deepcopy(self)
            if verbose:
                print "rep number " + str(i) + ' of ' + str(num_reps)
            res = jackknife_tree.evaluateSingleBalancedUPGMA(num_to_keep, \
                            Weight, norm)
            for key, ct in res.items():
                if key in result:
                    result[key] += ct
                else:
                    result[key] = 1
        return result

    def evaluateSingleBalancedUPGMA(self, num_to_keep, Weight, norm=False):
        """determine single iteration supported nodes in balanced env trees
        
        num_to_keep: The number of sequences to keep during each jackknife rep
        Weight: if False, makes tree using unweighted UniFrac; if True, uses
                weighted UniFrac
        norm: Whether or not to normalize the weighted UniFrac values by
                dividing by the average distance of seqs from the root"""
        result = {}
        self.balanceTreeEnvs(num_to_keep, Weight)
#hhere
        cluster = self.clusterEnvsUPGMA(Weight=Weight, norm=norm)
        cluster.nameUnnamedNodes()
        clusterArray, node_list = cluster.makeTreeArray()
        for i in clusterArray:
            key = tuple(i)
            if key in result:
                result[key] += 1
            else:
                result[key] = 1
        return result

    def assignJackknifeCounts(self, env_counts, pop_size=1.0):
        """
        Takes dictionary of environment:count and sets JackknifeCount attribute

        Sets count to 0 if no environment count available
        """
        for node in self.traverse():
            if node.Children:
                cur_key = tuple(node.Envs)
                if cur_key in env_counts:
                    node.JackknifeCount = env_counts[cur_key]
                    node.JackknifeFraction = \
                            (node.JackknifeCount * 1.0) / pop_size 
                    node.Data = None
                else:
                    node.JackknifeCount = 0
                    node.JackknifeFraction = 0.0
                    node.Data = None

    def parse_balanced_result(self, result, Weight=False):
        """parses output of evaluateBalancedUPGMA into tab-del text

        output_all lists all of the nodes that were found and the # of times
        ouput_nodes_str outputs info only for nodes in the original tree"""
        
        output_all = []
        output_nodes = {}
        output_nodes_str = []
        self.assignInternalEnvs()
        all_envs = [self.EnvsNames[i] for i in nonzero(self.Envs)]
        cluster = self.clusterEnvsUPGMA(Weight)
        cluster.nameUnnamedNodes()
        clusterArray, node_list = cluster.makeTreeArray()
        for i in clusterArray:
            key = tuple(i)
            if key in result:
                output_nodes[key] = result[key]
        for i in result:
            node = [all_envs[j] for j in nonzero(i)]
            output_all.append(str(node) + '\t' + str(result[i]))
            
        for i in output_nodes:
            node = [all_envs[j] for j in nonzero(i)]
            output_nodes_str.append(str(node) + '\t' + str(output_nodes[i]))
        return ('\n').join(output_all), ('\n').join(output_nodes_str)
        
#==============================================================
    
    def _find_sig_nodes_a(self, critical_p, fxn):
        """Returns nodes in self with a p-value below critical_p.

        fxn = the algorithm function used to test for significance
            EnvsNode.nodeChiSquare for Chi Square
            EnvsNode.nodeG_fit for G test of Significance
        
        Used with functions that:
            1) take the node as an input
            2) Can be calculated without removing the node from the tree
            3) Do not need pop-size as an input parameter

        critical_p: p_value under which nodes will be reported

        """
        node_list = [node for node in self.traverse() if node.Parent and \
                     node.Children]
        sig_nodes, all_nodes = self._find_sig_nodes_list_a\
                (node_list, fxn, critical_p) 
        return sig_nodes, all_nodes

    def _find_sig_nodes_list_a(self, nodes, fxn, critical_p):
        """Returns nodes in list, nodes, with a p-value below critical_p
        
        nodes: a list of nodes on which the test will be run
        
        fxn = the algorithm function used to test for significance
            EnvsNode.nodeChiSquare for Chi Square
            EnvsNode.nodeG_fit for G test of Significance

        critical_p: p-value under which nodes will be reported
        """
        sig_nodes = []
        all_nodes = []
        for node in nodes:
            fxn(self, node)
            #correct for multiple comparisons
            p_val = node.p * len(nodes)
            if p_val > 1.0:
                p_val = 1.0
            if p_val <= critical_p:
                sig_nodes.append((node, str(p_val)))
            all_nodes.append((node, str(p_val)))
        return sig_nodes, all_nodes

    def _find_sig_nodes_b(self, fxn, critical_p, pop_size, verbose=False):
        """Returns nodes in self with a p-value below critical_p
        
        Reports a list of nodes that have a significant p-value
                with only that node in the tree.

            fxn = the algorithm function used to test for significance
                EnvsNode.phylogeneticTestTree 
                EnvsNode.uniFracTree
                EnvsNode.uniFracWeightedTree
                
            :critical_p = p_value under which nodes will be reported
            :pop_size = the number of replication used to generate the p_value
                each time
        """
        nodes = [node for node in self.traverse() if node.Children]
        sig_nodes, all_nodes = self._find_sig_nodes_list_b(nodes, fxn, \
                critical_p, pop_size)
        return sig_nodes, all_nodes

    def _find_sig_nodes_list_b(self, nodes, fxn, critical_p, pop_size, \
            verbose=False):
        """Returns nodes in list, nodes with a p-value below critical_p
        
        Reports a list of nodes that have a significant p-value
                with only that node in the tree. Performs the 
                randomization for the whole tree but returns the p-value
                for the node alone

            fxn = the algorithm function used to test for significance
                EnvsNode.phylogeneticTestTree 
                EnvsNode.uniFracTree
                EnvsNode.uniFracTreeWeight
            
            :nodes: a list of nodes to evaluate
            :critical_p = p_value under which nodes will be reported
            :pop_size = the number of replications used to generate the p_value
                each time

        """
        sig_nodes = []
        all_nodes = []
        num_comp = len(nodes)
        for node in nodes:
            p_val, p_val_str = self._get_corrected_ps(node, fxn, pop_size, \
                        num_comp, verbose)
            if p_val <= critical_p:
                sig_nodes.append((node, p_val_str))
            all_nodes.append((node, p_val_str))
        return sig_nodes, all_nodes

    def _get_corrected_ps(self, node, fxn, pop_size, num_comp, verbose):
        """returns a p_val and a p_val string with p_val corrections

        called by _find_sig_nodes_list_b
        """
        fxn(self, pop_size, node)
        if node.p is None:
            return 1.0, "NA"
        if node.p == 0.0:
            p_val = (1/float(pop_size)) * num_comp
            if p_val > 1.0:
                p_val = 1.0
            p_val_str = '<=' + str(p_val)
        else:
            p_val = node.p * num_comp
            if p_val > 1.0:
                p_val = 1.0
            p_val_str = str(p_val)
        if verbose:
            print node.Data, p_val_str
        return p_val, p_val_str

    def sigNodeGrouper(self, sig_nodes):
        """determines which sig nodes are descendants of the others"""
        #dict for storing lists of sig nodes that are a dec of each sig node 
        node_groups = {}
        nodes = [node[0] for node in sig_nodes]
        for sig_node in nodes:
            # get all non-terminal descendants of sig_node
            decs = [node for node in sig_node.traverse()if node.Children]

            #list of sig nodes that are descendant of sig_node
            
            sig_decs = [n.Data for n in nodes if n in decs]

            node_groups[sig_node.Data] = sig_decs

        return node_groups
                    
    def generateWithSigOutput(self, sig_nodes, all_nodes, delimiter = '\t'):
        """generates text output for _find_sig_nodes_b"""
        output_str = ""
        #make label line
        label_line_lis = ["Node", "p-value", "num taxa"]
        label = delimiter.join(label_line_lis)

        #add results for significant nodes
        output_str += "Results for significant nodes\n"
        output_str += label
        output_str += '\n'
        for node in sig_nodes:
            data = [node[0].Data, node[1], str(len(node[0].TerminalDescendants))]
            s = delimiter.join(data)
            output_str += s
            output_str += '\n'
            
        #add node group information
        output_str += '\n\n'
        output_str = output_str + "Significant node groups\n" +\
            "Lists all significant nodes that are children of " +\
            "a significant node.\n\n"
        node_groups = self.sigNodeGrouper(sig_nodes)
        for group in node_groups:
            output_str += group
            output_str += '= \t'
            output_str += ', '.join(node_groups[group])
            output_str += '\n\n'

        #add results for all nodes
        output_str += "Results for all nodes\n"
        output_str += label
        output_str += '\n'
        for node in all_nodes:
            data = [node[0].Data, str(node[1]), \
                        str(len(node[0].TerminalDescendants))]
            s = delimiter.join(data)
            output_str += s
            output_str += '\n'

        #add Newick strings for all nodes
        output_str += '\n\nNewick string for all significant nodes:\n\n'
        for node in sig_nodes:
            output_str = output_str + node[0].Data + " = " + str(node[0]) + ';' +\
                         "\n\n"
        return output_str
            
#==============================================================

    def find_nodes_threshold(self, threshold):
        """returns nodes formed if tree was cut at a given distance from root

        the distance is defined by threshold"""
        #find nodes with a distance from root that is >= threshold but whose
        #parent is before the threshold
        result = []
        for node in self.traverse():
            if node is None or node.Parent is None:
                continue
            if self.distance(node) >= threshold and self.distance(node.Parent)\
                    < threshold:
                   result.append(node)
        return result
       
    def evaluate_threshold_nodes(self, threshold, fxn, critical_p=0.05, \
            pop_size=1000, min_leafs=4):
        """returns nodes at a given distance from root that have a sig. p-value

        threshold: a distance from the root at which to cut the tree. 
        fxn = the algorithm function used to test for significance
                EnvsNode.phylogeneticTestTree 
                EnvsNode.uniFracTree
                EnvsNode.uniFracWeightedTree
                EnvsNode.nodeChiSquare for Chi Square
                EnvsNode.nodeG_fit for G test of Significance
                
            
        :critical_p = p_value under which nodes will be reported
        :pop_size = the number of replication used to generate the p_value
                each time
        :min_leafs = min number of descendants that a node must have to be 
                evaluated
                
        sig_nodes is a list of (node, p-value) tuples for the nodes with 
            corrected p-values that are below the critical p
        all_nodes is a list of (node, p-value) tuples for all  of the evaluated
            nodes
            
        ***Currently only use for G-fit - PhylogeneticTestTree and UniFrac
        now randomize sequences among the whole tree when evaluating a single
        node and we are not sure if this is the best method. This functionality
        is also not well tested."""
        node_list = self.find_nodes_threshold(threshold)
        self.clearInternalEnvs
        self.assignInternalEnvs(Weight=True)
        new_list = [node for node in node_list if sum(node.Envs) \
                >= min_leafs]
        if fxn == EnvsNode.phylogeneticTestTree or fxn == EnvsNode.uniFracTree \
                or fxn == EnvsNode.uniFracWeightedTree:
            sig_nodes, all_nodes = self._find_sig_nodes_list_b\
                    (new_list, fxn, critical_p, pop_size)
        elif fxn == EnvsNode.nodeChiSquare or fxn == EnvsNode.nodeG_fit:
            sig_nodes, all_nodes = self._find_sig_nodes_list_a\
                    (new_list, fxn, critical_p)
        else:
            raise ValueError, "Invalid algorithm %s" % str(fxn)
        return sig_nodes, all_nodes

#===============================================================
#general use functions

    def mapNodeP(self):
        """maps the Data attribute of each node in self to its p attribute

        p is the p-value that can be assigned with a variety of methods"""
        p_val_info = {}
        for node in self.traverse():
            if node.Children:
                p_val_info[node.Data] = node.p
        return p_val_info
    
    def findCloseEnvs(self, node, threshold):
        """find environments assigned to close relatives of node

        node is the node to look for relatives of
        threshold is the branch length threshold to include as a relative
        """
        relatives = self.findCloseRelatives(node, threshold)
        related_envs = []
        if relatives:
            for node in relatives:
                env_list = [node.EnvsNames[i] for i in nonzero(node.Envs)]
                for i in env_list:
                    if i not in related_envs:
                        related_envs.append(i)
        return related_envs
                                
#=====================================================
#Future TreeNode Methods

    def reassignData(self, taxa_names):
        """reassigns the Data property of terminal descendants based on a dict:

        Takes dict taxa_names with new taxa names indexed by the old taxa names
        and reassigns the Data property of terminal descendants. This is for
        Nexus formatted tree files where each node is referenced by a number
        but also has a taxa name.
        """
        for term_dec in self.TerminalDescendants:
            if term_dec.Data in taxa_names:
                term_dec.Data = taxa_names[term_dec.Data]

    def removeNodeAndChildlessParents(self, node):
        """removes a specified node from the tree and parent if childless"""
        #get the parent for the node to remove
        parent = node.Parent

        #remove the node and its otherwise childless ancestors
        while node.Parent is not None: #while not at the root
            #dissociates node from tree by breaking link with parent
            node.Parent = None
            if parent: #if the parent has other children don't remove it
                break
            else:
                node = parent
                parent = node.Parent

    def removeNodeAndOnlyChildParents(self, node):
        """removes a specified node from the tree and parent if childless

        Also removes parent node if it only has one child and makes necessary
        adjustments to parent-child relationships and node branchlengths"""
        #get the parent for the node to remove
        parent = node.Parent
    
        #remove the node and its otherwise childless ancestors
        while node.Parent: #while not at the root
            parent_bl = parent.BranchLength
            #dissociates node from tree by breaking link with parent
            node.Parent = None
            #if parent has more than one additional children don't remove it
            if len(parent.Children) > 1: 
                break
            #if parent has one additional child remove it and make corrections
            if len(parent.Children) == 1 and parent.Parent:
                grandparent = parent.Parent
                #remove parent from tree
                parent.Parent = None
                #reattach remaining child to its grandparent
                child = parent.Children[0]
                child.Parent = grandparent
                #if there is branchlength information for the removed parent,
                #add it to the branchlength of the remaining child
                if parent_bl:
                    child.BranchLength += parent_bl
            #if parent has one child and grandparent is root,
            #attach child to root
            elif len(parent.Children) == 1 and not parent.Parent:
                child = parent.Children[0]
                grandchildren = child.Children
                child.Parent = None
                for grandchild in grandchildren:
                    grandchild.Parent = parent
                parent.BranchLength = child.BranchLength
            else:
                node = parent
                parent = node.Parent
    
    def keepTermDecsData(self, keep):
        """Keeps only the terminal descendants whose Data property are in keep

        keep is a list of names of the Terminal Descendants to keep"""

        for node in self.TerminalDescendants:
            if node.Data not in keep:
                self.removeNodeAndOnlyChildParents(node)
    
    def getSubtreeByLCA(self, taxa1, taxa2):
        """Splits the tree at the LCA of taxa1 and taxa2 and returns the LCA.

        This function makes it possible to calculate statistics
        a) without a clade of interest included or
        b) only including the clade of interest

        the function finds the Last Common Ancestor of two nodes that are
        identified by their Data property. By setting the Parent of the
        LCA to none, the input tree now does not include that clade.
        The LCA can be evaluated on its own as a new tree with the LCA
        as the root. """
        #identify the Node objects from their Data attributes
        for dec in self.TerminalDescendants:
            if dec.Data == taxa1:
                node = dec
            if dec.Data == taxa2:
                other = dec
        #find the LCA of the 2 nodes
        LCA = node.lastCommonAncestor(other)
        #remove clade from tree
        LCA.Parent = None
        return LCA

    def getSubtreeByData(self, node):
        """splits the tree at the node"""
        node.Parent = None
        return node
           
    def selectNodeByLeafCount(self, num_leaves):
        """returns a node with the specified number of terminal descendants"""
	for node in self.traverse():
            if len(node.TerminalDescendants) == num_leaves:
                if not node.Data:
                    node.Data = 'sel_node'
                return node
                break

    def mapNodeNames(self):
        """maps the Data attribute of each node obj in self to the node object
	"""
        result = {}
        for node in self.traverse():
            if node.Data:
                result[node.Data] = node
            else:
                pass
        return result

#==============================================
#Future PhyloNode methods

    def assignBranchLength(self, branch_lengths):
        """Sets BranchLength property of Nodes with info in branch_lengths dict

        branch_lengths is a dictionary that maps the name for each node to a
        tuple with the parent and the branchlength that is associated with it

        This method is used to populate Branch Length information from a 
        PAUP log file which reports  branch lengths
        in a table."""
        #post-order tree traversal
        for node in self.traverse(False, True):
            if node == self.Root:
                pass
            else:
                node.BranchLength = branch_lengths[node.Data][1]
                parent = node.Parent
                parent.Data = branch_lengths[node.Data][0]

    def setBranchLength(self, value=1.0):
        """Sets the BranchLength prop of all nodes except the root to value"""

        for node in self.traverse():
            if node.Parent is not None:
                node.BranchLength = value

    def totalBranchLength(self, env1=None, env2=None):
        """returns total branch length of the tree
        
        optionally takes a list of two environment names (envs) and only counts
        branches in those envs"""
        if not env1 and not env2:
            lengths = self.envBranchLengthsTree()
        else:
            totlengths = self.envBranchLengthsTree()
            lengths = self.envBranchLengthsPair(env1, env2, totlengths)
        totalbl = sum(lengths.values())
        return totalbl

    def makeNodalDistanceMatrix(self, default=None):
        """makes a dict2d of distances between Terminal Descendants"""
        matrix = []
        nodes = self.TerminalDescendants
        node_names = [node.Data for node in nodes]
        node_names.sort()
        for index, node in enumerate(nodes):
            node1 = node
            index2 = index + 1
            while index2 < len(nodes):
                node2 = nodes[index2]
                index2 += 1
                dist = node1.distance(node2)
                matrix.append((node1.Data, node2.Data, dist))
                matrix.append((node2.Data, node1.Data, dist))
        return Dict2D(matrix, RowOrder=node_names, ColOrder=node_names, \
                        Pad=True, Default=default)
                                
    def findCloseRelatives(self, node, threshold):
        """finds seqs in the tree within a branch length threshold

        node is the tree object to find the relatives of
        threshold is a branch length threshold
        """
        result = []
        #first get the deepest ancestor node that is still within threshold
        curr = node
        while node.distance(curr) <= threshold:
            if curr.Parent:
                curr = curr.Parent
            else:
                break
        #append only the children that are within threshold to result
        if curr.Children:
            for dec in curr.TerminalDescendants:
                if dec.distance(node) <= threshold:
                    result.append(dec)
        return result
                                                                                    

#===================================
#functions for determining sequence saturation

def get_sat_curve_data_envs(tree, num_runs, envs=None, stride=1):
    """for each env in envs, generate a saturation curve
    
    The taxa are removed from a particular environment leaving the rest
    of the tree intact. It can thus show the amount of diversity that is
    unique to a particular environment/ environment type
    
    envs is a list or set of environment names. If it is "all" does for
    all envs

    note: this does not yet behave correctly for trees in which the same taxa 
    may be assigned to multiple environments
    """
    
    result = []
    if not envs:
        envs = [tree.EnvsNames[i] for i in nonzero(tree.Envs)]
    for env in envs:
        env_result = {}
        num_taxa_env = len(tree.getEnvTermDecs(env))
        for r in range(num_runs):
            num_taxa = num_taxa_env
            new_tree = deepcopy(tree)
            num_rem = 0
            if str(num_rem) not in env_result:
                env_result[str(num_rem)] = []
            env_result[str(num_rem)].append(new_tree.totalBranchLength())
            while num_taxa >= stride:
                num_taxa -= stride
                num_rem += stride
                new_tree.removeEnvTermNodes(env, stride)
                key = str(num_rem)
                if key not in env_result:
                    env_result[key] = []
                env_result[key].append(new_tree.totalBranchLength())
        result.append(env_result)
    return result, envs

def get_sat_curve_data_tree(tree, num_runs, stride=1, Verbose=False):
    """remove any sequence from a tree sequentially and reports total bl.

    Used to generate a saturation curve for a whole tree without regard to
    environment assignments"""
    result = {}
    num_taxa_tree = len(tree.TerminalDescendants)
    for r in range(num_runs):
        if Verbose:
            print "run number " + str(r) + 'of' + str(num_runs)
        num_taxa = num_taxa_tree
        new_tree = deepcopy(tree)
        if str(num_taxa) not in result:
            result[str(num_taxa)] = []
        result[str(num_taxa)].append(new_tree.totalBranchLength())
        while num_taxa > stride:
            num_taxa -= stride
            new_tree.removeEnvTermNodes(None, stride)
            if str(num_taxa) not in result:
                result[str(num_taxa)] = []
            result[str(num_taxa)].append(new_tree.totalBranchLength())
    return result

def generate_sat_curve_output(sat_curve_info, delimiter = '\t', \
        num_taxa=None, complete=True):
    """turns output of get_sat_curve_data into a string for output"""
    output = []
    #make label line
    label_line_lis = ["num_taxa", "average_bl", delimiter, "raw_data"]
    s = delimiter.join(label_line_lis)
    output.append(s)

    #add data lines
    #get the number of taxa in the tree if it was not given as input
    if not num_taxa:
        num_taxa = max([int(key) for key in sat_curve_info.keys()])
    #for each possible number of taxa, make a line, if no calculations
    #were made the line will be left blank. This is so that different
    #datasets can be concatenated and compared
    for i in range(num_taxa + 1):
        line = str(i) + delimiter
        if str(i) in sat_curve_info:
            average = sum(sat_curve_info[str(i)])/len(sat_curve_info[str(i)])
            line = line + str(average) + delimiter + delimiter + delimiter
            #only add the raw data (in addition to the average, if the bl
            #is true
            if complete:
                add = [str(val) for val in sat_curve_info[str(i)]]
                line += (delimiter).join(add)
        output.append(line)

    return output

#==========================================================
#Generating tree objects/ parsing files

def parse_newick(tree_f):
    """extracts the dnd from a Newick from ARB tree file
    """
    dnd = {}
    in_tree = False
    dnd["tree_1"] = []
    for line in tree_f:
            if line.startswith("("):
                in_tree = True
            if in_tree:
                dnd["tree_1"].append(line)
            if line.endswith(";"):
                in_tree = False
    if dnd["tree_1"]:
        return dnd
    else:
        raise RecordError, "not a valid Newick file."

def make_env_dict_text_Weight(Trans_table):
    """makes an env dict from tab delimited text that includes abundance info
    
    This option should be used for a weighted analysis. The abundance is 
    information of how many of a particular sequence was found as determined
    with RFLP or based on sequence information
    
    format: taxa_name \t environment \t abundance

    returns a dictionary mapping taxa to an array listing whether the taxa is
    present in the different environments represented in the tree

    also returns a list of the total environments represented in the tree
    that cooresponds to the arrays in envs_dict
    """

    #get a list of all of the envs in the file
    all_envs = []
    line_info = []
    env_dict = {}

    for line in Trans_table:
        info = line.split('\t')
        if len(info) == 3:
            taxa_name = info[0].strip()
            env_name = info[1].strip()
            weight = float(info[2].strip())
            line_info.append((taxa_name, env_name, weight))
            if env_name not in all_envs:
                all_envs.append(env_name)
    all_envs.sort()
    for item in line_info:
        taxa_name = item[0]
        env_name = item[1]
        weight = item[2]
        if taxa_name not in env_dict:
            #set the taxa_name to an array with 0 for every env in the tree
            env_dict[taxa_name] = zeros(len(all_envs), Float)
            #set the index for the environment to the abundance
            env_dict[taxa_name][all_envs.index(env_name)] = weight
        else:
            #set the index for the environment to the abundance
            env_dict[taxa_name][all_envs.index(env_name)] = weight
    all_envs.sort()
    return env_dict, all_envs
    
def make_env_dict_text(Trans_table):
    """makes an env dict from tab delimited text

    format = taxa_name \t environment \t environment \n

    returns a dictionary mapping taxa to an array listing whether the taxa is
    present in the different environments represented in the tree

    also returns a list of the total environments represented in the tree
    that cooresponds to the arrays in envs_dict
    """
    #make an env_dict mapping taxa name to a set of environment information
    env_dict = {}
    for line in Trans_table:
        info = line.split('\t')
        if len(info) >= 2:
            name = info.pop(0).strip()
            env_dict[name] = Set([])
            for i in info:
                env_dict[name] = env_dict[name] | Set([i.strip()])
        else:
            pass
    #get the union of all of the environment sets to get the complete set
    #of environments in the tree
    all_envs = Set([])
    for i in env_dict:
        all_envs = all_envs.union(env_dict[i])
    #change to a sorted list and use to set EnvsNames
    all_envs = list(all_envs)
    all_envs.sort()

    #replace the sets in envs_dict  with arrays representing the envs  
    for i in env_dict:
        #initialize an envs array
        node_envs = array(zeros((len(all_envs),),Float))
        #populate node_envs based on env set (i)
        for index, tree_env in enumerate(all_envs):
            if tree_env in env_dict[i]:
                node_envs[index] = 1
        env_dict[i] = node_envs
    return env_dict, all_envs

def convert_envDict_to_tabDel(env_dict):
    """converts an env_dict to a tab-delimited string for import"""
    result = []
    for i in env_dict:
        line = i + '\t'
        for env in env_dict[i]:
            line = line + env + '\t'
        line += '\n'
        result.append(line)
    return result

def balanced_tree(levels, Constructor):
    """creates a strictly bifurcating tree at specified # of levels
    
    Constructor is the type of Tree object to use
        TreeNode
        PhyloNode
        EnvsNode
        EnveNodeICR"""
    root = Constructor()
    for i in range(levels):
        for t in root.traverse():
            if not t:
                t.append(Constructor())
		t.append(Constructor())
    return root

def make_tree_list(tree_lines, env_lines, env_format):
    """makes a list of tree objects from open file objects

    tree_lines: contains one tree in Newick format per line. 
    env_lines: open file object for tab delimited envs file
    envs_format: format of input file:
        TabDel
        TabDelWeight
"""
    if env_format == 'TabDel':
        envs, all_envs = make_env_dict_text(env_lines)
    elif env_format == 'TabDelWeight':
        envs, all_envs = make_env_dict_text_Weight(env_lines)
    else:
        return ValueError, "Not a valid format code."
    trees = []
    dnd = []
    def is_tail_line(line):
            return line.endswith(';')
    NewickFinder = TailedRecordFinder(is_tail_line)
    for dnd in NewickFinder(tree_lines):
        tree = DndParser(dnd, constructor = EnvsNode)
        tree.assignTerminalEnvs(envs, all_envs)
        trees.append(tree)
    return trees

def batch_pairwise_unifrac_values(trees, Weight=False, norm=False):
    """computes pairwise UniFrac values for all tree objects in trees

    trees: a list of tree objects
    Weight: if False performs un-weighted UniFrac, if True performs weighted
    
    for all trees performs pairwise UniFrac comparisons, returns a dict
    of dicts with a list of values for each env combo"""
    
    result = {} 
    for tree in trees:
        tree_result = tree.makeEnvMatrix(Weight=Weight, norm=norm)
        for env1 in tree_result:
            if env1 not in result:
                result[env1] = {}
            for env2 in tree_result[env1]:
                if env2 not in result[env1]:
                    result[env1][env2] = [tree_result[env1][env2]]
                else:
                    result[env1][env2].append(tree_result[env1][env2])
    return Dict2D(result)
        
def average_pairwise_unifrac_values(trees, default=None, Weight=False, norm=False):
        """makes a dict2D with average totalUnique for pairwise env combos for tree set

        trees: a list of tree objects

        default: Pads the dict2D object with this value

        Weight: if False, does not account for abundance info; if true it does

        norm: if True, normalizes Weighted UniFrac values so they are a value
        between 0 and 1
        """
        average_values = []
        batch_values = batch_pairwise_unifrac_values(trees, Weight, norm)
        batch_values.setDiag(None)
        for env1 in batch_values:
                for env2 in batch_values[env1]:
                        if batch_values[env1][env2]:
                                all_values = batch_values[env1][env2]
                                average = sum(all_values)/len(all_values)
                                average_values.append((env1, env2, average))
        envs = [trees[0].EnvsNames[i] for i in nonzero(trees[0].Envs)]
        return Dict2D(average_values, RowOrder = envs, ColOrder=envs, Pad=True, \
                        Default=default)
                                               
                                
def uniFracP_uncertainty(trees, Weight=False, env1=None, env2=None):
    """uniFracP_uncertainty takes a list of trees as input and returns p-value

    trees: a list of tree objects
    Weight: if True - a weighted analyis is done

    This function takes a list of tree objects as input and returns a single
    p-value. The list of tree objects represent equally likely trees returned
    from a liklihood analysis or parsimony or a group of bootstrapped trees.
    For each tree, the real UniFrac value is first calculated and then a 
    UniFrac value is calculated after randomizing the environment assignments.
    The two lists are then evaluated with a paired T-test to determine if
    the distributions of real and random values are significantly different
    from each other. A one-tailed test is used since you would expect to 
    have _more_ Unique branchlength if the environments are different.
    """
    real_vals = []
    random_vals = []
    for tree in trees:
        if not Weight:
            random_tree = deepcopy(tree)
            real_vals.append(tree.totalUnique())
            random_tree._make_randomized()
            random_vals.append(random_tree.totalUnique())
    t, prob = t_paired(real_vals, random_vals, 'high')
    return t, prob

def batch_cluster_envs_UPGMA(trees, Weight=False, norm=False):
    """returns a tree clustering the environments with UPGMA based on a tree set.
        
        trees: a list of EnvsNode objects
        
        Weight: Determines whether to use normal(False) or Weighted(True)
                UniFrac
               
        norm: Determine whether to normalize the raw weighted UniFrac values
                (True) or not (False)"""

    #make a distance matrix for envs and pad with a large number
    env_matrix = average_pairwise_unifrac_values(trees, BIG_NUM, Weight, norm)
    #turn dict2D into Numeric array and make list of phylonode objects
    matrix_a, node_order = inputs_from_dict2D(env_matrix)
    return UPGMA_cluster(matrix_a, node_order,BIG_NUM)

def batch_cluster_envs_PCA(trees, Weight=False, norm=False):
    """Principal coordinates analysis of envs using average UniFrac dist. for a tree set

        returns tab delimited text that can be analyzed in Excel        
        """
        
    dis_mat = average_pairwise_unifrac_values(trees, 0.0, Weight, norm)
    matrix, node_order = inputs_from_dict2D(dis_mat)
    names = [node.Data for node in node_order]
    point_matrix, eigvals = principal_coordinates_analysis(matrix)
    return output_pca(point_matrix, eigvals, names)


#def uniFracP_uncertainty_pairwise(trees, Weight=False, env1=None, env2=None):
#    """uniFracP_uncertainty takes a list of trees as input and returns p-value
#
#    trees: a list of tree objects
#    Weight: if True - a weighted analyis is done
#
#    This function takes a list of tree objects as input and returns a single
#    p-value. The list of tree objects represent equally likely trees returned
#    from a liklihood analysis or parsimony or a group of bootstrapped trees.
#    For each tree, the real UniFrac value is first calculated and then a 
#    UniFrac value is calculated after randomizing the environment assignments.
#    The two lists are then evaluated with a paired T-test to determine if
#    the distributions of real and random values are significantly different
#    from each other. A one-tailed test is used since you would expect to 
#    have _more_ Unique branchlength if the environments are different.
#    """
#
#    matrix = []
#    #get the number of comparisons being made
#    example_tree = trees[0]
#    num_envs = len([example_tree.EnvsNames[i] for i in \
#            nonzero(example_tree.Envs)])
#    num_comp = sum([i for i in range(1,num_envs)])
#    #get p-values for all pair combinations
#    for env1, env2 in self.getEnvPairs():
#        real_vals = []
#        random_vals = []
#        for tree in trees: 
#            if not Weight:
#                random_tree = deepcopy(tree)
#                real_vals.append(tree.totalUnique())
#                random_tree._make_randomized()
#                random_vals.append(random_tree.totalUnique())
#        t, prob = t_paired(real_vals, random_vals, 'high')
#    return t, prob

    
